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Abstract 

We consider the probabilistic numerical scheme for fully nonlinear PDEs suggested 
in |12j , and show that it can be introduced naturally as a combination of Monte Carlo 
and finite differences scheme without appealing to the theory of backward stochastic dif- 
ferential equations. Our first main result provides the convergence of the discrete-time 
approximation and derives a bound on the discretization error in terms of the time step. 
An explicit implementable scheme requires to approximate the conditional expectation 
operators involved in the discretization. This induces a further Monte Carlo error. Our 
second main result is to prove the convergence of the latter approximation scheme, 
and to derive an upper bound on the approximation error. Numerical experiments are 
performed for the approximation of the solution of the mean curvature flow equation in 
dimensions two and three, and for two and five-dimensional (plus time) fully-nonlinear 
Hamilton- Jacobi-Bellman equations arising in the theory of portfolio optimization in 
financial mathematics. 
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1 Introduction 

We consider the probabilistic numerical scheme for the approximation of the solution of a 
fully-nonlinear parabolic Cauchy problem suggested in [12]. In the latter paper, a repre- 
sentation of the solution of the PDE is derived in terms of the newly introduced notion of 
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second order backward stochastic differential equations, assuming that the fuhy-nonhnear 
parabohc Cauchy problem has a smooth solution. Then, similarly to the case of backward 
stochastic differential equations which are connected to semi-linear PDEs, this representa- 
tion suggests a backward probabilistic numerical scheme. 

The representation result of can be viewed as an extension of the Feynman-Kac rep- 
resentation result, for the linear case, which is widely used in order to approach the numer- 
ical approximation problem from the probabilistic viewpoint, and to take advantage of the 
high dimensional properties of Monte Carlo methods. Previously, the theory of backward 
stochastic differential equations provided an extension of these approximation methods to 
the semi- linear case. See for instance Chevance |13j . El Karoui, Peng and Quenez |18| . Bally 
and Pages [;2J, Bouchard and Touzi [9j and Zhang [32]. In particular, the latter papers pro- 
vide the convergence of the "natural" discrete-time approximation of the value function 
and its partial space gradient with the same L? error of order ^/h, where h is the length of 
time step. The discretization involves the computation of conditional expectations, which 
need to be further approximated in order to result into an implementable scheme. We refer 
to [2], [9] and [20] for an complete asymptotic analysis of the approximation, including the 
regression error. 

In this paper, we observe that the backward probabilistic scheme of [12] can be introduced 
naturally without appealing to the notion of backward stochastic differential equation. This 
is shown is Section [2] where the scheme is decomposed into three steps: 

(i) The Monte Carlo step consists in isolating the linear generator of some underlying 
diffusion process, so as to split the PDE into this linear part and a remaining nonlinear 
one. 

(ii) Evaluating the PDE along the underlying diffusion process, we obtain a natural discrete- 
time approximation by using finite differences approximation in the remaining nonlinear 
part of the equation. 

(iii) Finally, the backward discrete-time approximation obtained by the above steps (i)-(ii) 
involves the conditional expectation operator which is not computable in explicit form. An 
implementable probabilistic numerical scheme therefore requires to replace such conditional 
expectations by a convenient approximation, and induces a further Monte Carlo type of 
error. 

In the present paper, we do not require the fully nonlinear PDE to have a smooth solution, 
and we only assume that it satisfies a comparison result in the sense of viscosity solutions. 
Our main objective is to establish the convergence of this approximation towards the unique 
viscosity solution of the fully- nonlinear PDE, and to provide an asymptotic analysis of the 
approximation error. 

Our main results are the following. We first prove the convergence of the discrete-time 
approximation for general nonlinear PDEs, and we provide bounds on the corresponding 
approximation error for a class of Hamilton-Jacobi-Bellman PDEs. Then, we consider the 
implementable scheme involving the Monte Carlo error, and we similarly prove a conver- 
gence result for general nonlinear PDEs, and we provide bounds on the error of approxi- 
mation for Hamilton-Jacobi-Bellman PDEs. We observe that our convergence results place 
some restrictions on the choice of the diffusion of the underlying diffusion process. First, a 
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uniform ellipticity condition is needed; we believe that this technical condition can be re- 
laxed in some future work. More importantly, the diffusion coefficient is needed to dominate 
the partial gradient of the remaining nonlinearity with respect to its Hessian component. 
Although we have no theoretical result that this condition is necessary, our numerical ex- 
periments show that the violation of this condition leads to a serious mis-performance of 
the method, see Figure [5j 

Our proofs rely on the monotonic scheme method developed by Barles and Souganidis 
[7] in the theory of viscosity solutions, and the recent method of shaking coefficients of 
Krylov [21], [25] and [26] and Barles and Jakobsen |6], [5] and [4|. The use of the latter 
type of methods in the context of a stochastic scheme seems to be new. Notice however, 
that our results are of a different nature than the classical error analysis results in the 
theory of backward stochastic differential equations, as we only study the convergence of 
the approximation of the value function, and no information is available for its gradient or 
Hessian with respect to the space variable. 

The following are two related numerical methods based on finite differences in the context 
of Hamilton- Jacobi-Bellman nonlinear PDEs: 

• Bonnans and Zidani [8] introduced a finite difference scheme which satisfies the crucial 
monotonicity condition of Barles and Souganidis [7] so as to ensure its convergence. 
Their main idea is to discretize both time and space, approximate the underlying 
controlled forward diffusion for each fixed control by a controlled local Markov chain 
on the grid, approximate the derivatives in certain directions which are found by 
solving some further optimization problem, and optimize over the control. Beyond 
the curse of dimensionality problem which is encountered by finite differences schemes, 
we believe that our method is much simpler as the monotonicity is satisfied without 
any need to treat separately the linear structures for each fixed control, and without 
any further investigation of some direction of discretization for the finite differences. 

• An alternative finite-differences scheme is the semi-Lagrangian method which solves 
the monotonicity requirement by absorbing the dynamics of the underlying state in 
the finite difference approximation, see e.g. Debrabant and Jakobsen |15j . Camilli 
and Jacobsen [11], Camilli and Falcone |10j . and Munos and Zidani |30j . Loosely 
speaking, this methods is close in spirit to ours, and corresponds to freezing the 
Brownian motion Wh, over each time step h, to its average order \/Ti. However it 
does not involve any simulation technique, and requires the interpolation of the value 
function at each time step. Thus it is also subject to the curse of dimensionality 
problems. 

We finally observe a connection with the recent work of Kohn and Serfaty [23j who provide 
a deterministic game theoretic interpretation for fully nonlinear parabolic problems. The 
game is time limited and consists of two players. At each time step, one tries to maximize her 
gain and the other to minimize it by imposing a penalty term to her gain. The nonlinearity 
of the fully nonlinear PDE appears in the penalty. Also, although the nonlinear penalty does 
not need to be elliptic, a parabolic nonlinearity appears in the limiting PDE. This approach 
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is very similar to the representation of [12] where such a parabohc envelope appears in the 
PDE, and where the Brownian motion plays the role of Nature playing against the player. 

The paper is organized as follows. In Section [2j we provide a natural presentation of 
the scheme without appealing to the theory of backward stochastic differential equations. 
Section [3] is dedicated to the asymptotic analysis of the discrete-time approximation, and 
contains our first main convergence result and the corresponding error estimate. In Section 
|4| we introduce the implementable backward scheme, and we further investigate the induced 
Monte Carlo error. We again prove convergence and we provide bounds on the approxi- 
mation error. Finally, Section [5] contains some numerical results for the mean curvature 
flow equation on the plane and space, and for a five- dimensional Hamilton-Jacobi-Bellman 
equation arising in the problem of portfolio optimization in financial mathematics. 

Notations For scalars a,b £ M, we write a Ab := min{a,6}, ay h := max{a, 6}, := 
max{— a, 0}, and a+ := max{a, 0}. By M(n, d), we denote the collection of all n x d matrices 
with real entries. The collection of all symmetric matrices of size d is denoted Sd, and its 
subset of nonnegative symmetric matrices is denoted by S^. For a matrix A G M(n, d), we 
denote by its transpose. For A,B£ M(n, d), we denote A-B := TvlA'^B]. In particular, 
for d = 1, A and B are vectors of M" and A ■ B reduces to the Euclidean scalar product. 

For a function u from [0, T] x M'^ to M, we say that u has g— polynomial growth (resp. 
a— exponential growth) if 



/ -a\x\\ ij. M \ 

sup j— j- < oo, (resp. sup e ' ' x)| < oo). 

For a suitably smooth function ip on Qt := (0, T\ X M^, we define 

II I / M 111 II ^) ~ ^01 

l^^loo := sup and |(y9|i := l^jjoo + sup ^. 

{i,x)6QT QtxQt {x - x') + \t- t'\ 2 

Finally, we denote the L^'— norm of a r.v. R by ||-R||p := (E[|i?|*'])"^''^. 

2 Discretization 

Let IX and a be two maps from x to and M((i, d), respectively. With a := aa^ . 
We define the linear operator: 

dt 2 

Given a map 

F : {t,x,r,p,'y) eR+ xR'^ xRxW^ xSd i — > F{x,r,p,-f) e R 

we consider the Cauchy problem: 

~C^v - F {■,v,Dv,D'^v) =0, on[0,T)xM^, (2.1) 
^(r,-) = g, on G R'^. (2.2) 
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Under some conditions, a stochastic representation for the solution of this problem was pro- 
vided in [12] by means of the newly introduced notion of second order backward stochastic 
differential equations. As an important implication, such a stochastic representation sug- 
gests a probabilistic numerical scheme for the above Cauchy problem. 

The chief goal of this section is to obtain the probabilistic numerical scheme suggested 
in [12] by a direct manipulation of (2.1 )-(2.2) without appealing to the notion of backward 
stochastic differential equations. 

To do this, we consider an M'^-valued Brownian motion on a filtered probability space 
($7, F, P), where the filtration F = {J-t,t G [0,7"]} satisfies the usual completeness condi- 
tions, and J^Q is trivial. 

For a positive integer n, let h := T/n, ti = ih, i = 0, . . . ,n, and consider the one step 
ahead Euler discretization 



X*'^' := X + fi{t, x)h + cj{t, x){Wt+h - Wt), 



(2.3) 



of the diffusion X corresponding to the linear operator . Our analysis does not require 
any existence and uniqueness result for the underlying diffusion X. However, the subsequent 
formal discussion assumes it in order to provides a natural justification of our numerical 
scheme. 



Assuming that the PDE (2.1) has a classical solution, it follows from Ito's formula that 



C^v{t,Xt)dt 



U 



where we ignored the difficulties related to local martingale part, and Kt-^c 
denotes the expectation operator conditional on {Xt- = x}. Since v solves the PDE (2.1), 
this provides 



v{ti,x) 



F{-,v,Dv,D'^v){t,Xt)dt 



By approximating the Riemann integral, and replacing the process X by its Euler dis- 
cretization, this suggest the following approximation of the value function v 



v^iT,.):=g and v^iU, x) := Th[v'']iti, x), 



where we denoted for a function : M_ 



E 



— > M with exponential growth: 

+ )] +hF{;Vhi^) {t,x), 



P^V'Ci,^) :=E[Z?V(t + /i,^f )], A; = 0,1, 2, V^i^ := {V^^i;,Vli;,Vli;) 



(2.4) 

(2.5) 
(2.6) 



and D'^ is the A;— th order partial differential operator with respect to the space variable x. 
The differentiations in the above scheme are to be understood in the sense of distributions. 
This algorithm is well-defined whenever g has exponential growth and -F is a Lipschitz 
map. To see this, observe that any function with exponential growth has weak gradient 
and Hessian because the Gaussian kernel is a Schwartz function, and the exponential growth 
is inherited at each time step from the Lipschitz property of F. 
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At this stage, the above backward algorithm presents the serious drawback of involving 
the gradient Dv^{ti^i, .) and the Hessian D'^v^{ti^i, .) in order to compute v^{ti,.). The 
following result avoids this difficulty by an easy integration by parts argument. 



Lemma 2.1. For every function ip : Qj- 



with exponential growth, we have: 



Vh^(ti,x) 
where Hh = {HI^',H[\H^Y and 



E 



ip{ti+i,Xj^''-^)Hh{ti,x 



(2.7) 



Proof. The main ingredient is the following easy observation. Let G be a one dimensional 
Gaussian random variable with unit variance. Then, for any function / : M — > M with 
exponential growth, we have: 



E[f{G)H\G)] = E[/('=)(G)], 



where /^-^^ is the A;— th order derivative of / in the sense of distributions, and is the 
one-dimensional Hermite polynomial of degree k. 

1 Now, let (p : M*^ — > M be a function with exponential growth. Then, by direct 



conditioning, it follows from (2.8) that 

E 

and therefore: 



d 



— (X^' )aji{t,x) 



E 



2 For i ^ j, it follows from (2.8) that 

d 



k=l 

d 



a{t,x)'E V^{X]l 



kl=l 



dxkdxi 



{Xff)cJij{t,x)aki{t,x) 



and for j = i: 



E 



This provides: 



^(X*''')cJ;i(t, x)aki{t, x) 

dxkOxi 



E 



^iX'r^)H^it,x) 



□ 



In view of Lemma 2.1 the iteration which computes v^{ti, .) out of .) in (2.4)-(2.5) 

does not involve the gradient and the Hessian of the latter function. 
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Remark 2.2. Clearly, one can proceed to different choices for the integration by parts in 
One such possibility leads to the representation of V ^s: 



Lemma 



2.1 



E 



ih/2) {h/2f 



This representation shows that the backward scheme (2.4) is very similar to the probabilistic 
numerical algorithm suggested in |12j . 



Observe that the choice of the drift and the diffusion coefficients ^ and a in the nonlinear 



PDE (2.1) is arbitrary. So far, it has been only used in order to define the underlying 



diffusion X. Our convergence result will however place some restrictions on the choice of 



the diffusion coefficient, see Remark 3.4 

Once the linear operator is chosen in the nonlinear PDE, the above algorithm handles 
the remaining nonlinearity by the classical finite differences approximation. This connection 



with finite differences is motivated by the following formal interpretation of Lemma 2.1 
where for ease of presentation, we set d = 1, /i = 0, and (t{x) = 1: 

• Consider the binomial random walk approximation of the Brownian motion Wt^ := 
kh, /c > 1, where {wj^j > 1} are independent random variables 

^ + 5_^^. Then, this induces the following approximation: 

il){t^ X + Vh) — '4){t, X — Vh) 



-.iWj, tk 



Vli^{t,x) :=E 



ij{t + h,x]f)H^ 



which is the centered finite differences approximation of the gradient. 

Similarly, consider the trinomial random walk approximation Wt^. := X^j=i Wj, '■= 
kh, k > 1, where {wj,j > 1} are independent random variables distributed as 

I {hVsh} + 4<5{0} + '^^{-v^}) > so that E[wJ 
this induces the following approximation: 



Vl^{t,x) 



E 



tp{t + h,X 



h ' 



= E[VF^] for all integers n < 4. Then, 
ij{t, x + VSh) - 2ip{t, x) + ij{t, X - V3h) 



3h 



which is the centered finite differences approximation of the Hessian. 

In view of the above interpretation, the numerical scheme studied in this paper can be 
viewed as a mixed Monte Carlo-Finite Differences algorithm. The Monte Carlo component 
of the scheme consists in the choice of an underlying diffusion process X. The finite dif- 
ferences component of the scheme consists in approximating the remaining nonlinearity by 



means of the integration- by-parts formula of Lemma 2.1 



3 Asymptotics of the discrete-time approximation 
3.1 The main results 

Our first main convergence results follow the general methodology of Barles and Souganidis 



[7j, and requires that the nonlinear PDE (2.1) satisfies a comparison result in the sense of 
viscosity solutions. 
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We recall that an upper-semicontinuous (resp. lower semicontinuous) function v (resp. 
v) on [0,T] X M'^, is called a viscosity subsolution (resp. supersolution) of (2.1) if for any 
{t, x) € [0, T) X M"' and any smooth function 93 satisfying 



= (v — (f)(t,x) = max (w — 99) resp. = (v — (f)(t,x) = min (v — ip) 
[0,r]xRrf V [0,T]xR'* 

we have: 

-C^ip- F{t,x,Vip{t,x)) < (resp. >) 0. 



Definition 3.1. We say that (2.1) has comparison for bounded functions if for any bounded 
upper semicontinuous subsolution v and any bounded lower semicontinuous supersolution v 
on [0, T) X M'', satisfying 

v{T,-)<v{T,-), 

we have v<v. 

Remark 3.2. Barles and Souganidis [T| use a stronger notion of comparison by accounting 
for the final condition, thus allowing for a possible boundary layer. In their context, a 
supersolution v and a subsolution v satisfy: 

mm{-C^v{T,x) - F{T,x,Vv{T,x)),v{T,x) - g{x)} < (3.1) 
max{-C^v{T,x) - F{T,x,Vv{T,x)),v{T,x) - g{x)} > 0. (3.2) 



We observe that, by the nature of our equation, (3.1) and (3.2) imply that the subsolution 
V < g and the supersolution v > g, i.e. the final condition holds in the usual sense, 
and no boundary layer can occur. To see this, without loss of generality we suppose that 



F{t,x,r,p,j) is decreasing with respect to r (see Remark 3.13). Let ip he a function 
satisfying 

= {v — (p)(T,x) = max (v — ip). 

[0,T]xR<^ 

Then define (fxit, •) = fit, •) + K{T — t) for K > 0. Then v — ipx also has a maximum at 



(T, x), and the subsolution property (3.1) implies that 

mm{-C^Lp{T,x)-F{T,x,Vip{T,x)) + K,v(T,x)-g{x)] < 0. 
For a sufficiently large K, this provides the required inequality v{T, x)—g{x) < 0. A similar 



argument shows that (3.1) implies that v — g > 0. 

In the sequel, we denote by i^^j Fp and F^ the partial gradients of F with respect to 
r, p and 7, respectively. We also denote by F~ the pseudo-inverse of the non-negative 
symmetric matrix i^^. We recall that any Lipschitz function is differentiable a.e. 

Assumption F (i) The nonlinearity F is Lipschitz- continuous with respect to (x,r, p, 7) 
uniformly in t, and \F{-, •, 0, 0, 0)|oo < 00; 

(ii) F is elliptic and dominated by the diffusion of the linear operator , i.e. 

VyF < a on M"* x M x M'' x S^; (3.3) 

(iii) Fp G Image(F^) and \F'jF~Fp\^ < +00. 
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Remark 3.3. Assumption F (iii) is equivalent to 



\^p\oo 



< oo where mp '■= min {Fp-w + w'^F^w} . (3.4) 



This is immediately seen by recalling that, by the symmetric feature of F^, any w G M'' has 
an orthogonal decomposition w = W1 + W2 G Ker(F^)©Image(-F^), and by the nonnegativity 



of F^: 



Fp ■ w + w F^w = Fp ■ wi + Fp ■ W2 + W2 F^W2 



\f^F-Fp + Fp.wi + \\{F^f'^ ■ Fp - F^^I^W2\\ 



Remark 3.4. The above Condition (3.3) places some restrictions on the choice of the 
linear operator in the nonlinear PDE (2.1 ). First, F is required to be uniformly elliptic, 



implying an upper bound on the choice of the diffusion matrix a. Since aa"^ G iSj", this 
implies in particular that our main results do not apply to general degenerate nonlinear 
parabolic PDEs. Second, the diffusion of the linear operator a is required to dominate the 
nonlinearity F which places implicitly a lower bound on the choice of the diffusion a. 



dv 



Example 3.5. Let us consider the nonlinear PDE in the one-dimensional case — ^ 



1 



2 



a V' 



lP'v~^ where < 6 < a are given constants. Then if we restrict the choice of the 
diffusion to be constant, it follows from Condition F that < o"^ < 6^, which implies that 
a?' < 3b^. If the parameters a and b do not satisfy the latter condition, then the diffusion 
a has to be chosen to be state and time dependent. 

Theorem 3.6 (Convergence). Let Assumption F hold true, and \^\i, < 00 and a is 



invertible. Also assume that the fully nonlinear PDE (2.1) has comparison for bounded 
functions. Then for every bounded Lipschitz function g, there exists a bounded function v 
so that 

— > V locally uniformly. 



In addition, v is the unique bounded viscosity solution of problem (2.1) -(2.2). 



Remark 3.7. Under the boundedness condition on the coefficients fi and a, the restriction 
to a bounded terminal data g in the above Theorem |3.6| can be relaxed by an immediate 
change of variable. Let 5 be a function with a— exponential growth for some a > 0. Fix 
some M > 0, and let p be an arbitrary smooth positive function with: 

p{x) = e"l^l for \x\ > M, 

so that both p{x)~^V p{x) and p{x)~^'S/'^ p{x) are bounded. Let 



u 



{t,x) := p{x)~'^v{t,x) for {t,x) e [0,T] x 



Then, the nonlinear PDE problem (2.1)-(2.2) satisfied by v converts into the following 
nonlinear PDE for u: 

- C^u- F {■,u,Du,D'^u) =0 on[0,r)xM'^ (3.5) 
v{T, ■)=g:= p-^g on 
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where 



F{t,x,r,p,^) := r^{x) ■ p p + -Tr [a{x) [r p ^V'^p + 2pp p^)]^ 



+p F it,x,rp,rV p + pp,rV p + 2pVp + p7j . 



Recall that the coefficients p and a are assumed to be bounded. Then, it is easy to see 
that F satisfies the same conditions as F. Since g is bounded, the convergence Theorem 
3.6 applies to the nonlinear PDE (3.5). □ 



Remark 3.8. Theorem 3.6 states that the inequality (3.3) (i.e. diffusion must dominate 
the nonlinearity in 7) is sufficient for the convergence of the Monte Carlo-Finite Differences 
scheme. We do not know whether this condition is necessary: 



• Subsection 3.4 suggests that this condition is not sharp in the simple linear case, 

• however, our numerical experiments of Section [5] reveal that the method may have a poor 
performance in the absence of this condition, see Figure [5j 

We next provide bounds on the rate of convergence of the Monte Carlo-Finite Differences 
scheme in the context of nonlinear PDEs of the Hamilton- Jacobi-Bellman type in the same 
context as [B]. The following assumptions are stronger than Assumption F and imply that 



the nonlinear PDE (2.1) satisfies a comparison result for bounded functions. 



Assumption HJB The nonlinearity F satisfies Assumption F (ii)-(iii), and is of the 
Hamilton-, Jacobi-Bellman type: 



-a • 7 + 6 • p + F{t, X, r,p, 7) 



inf {£"(t,rc,r,p,7)} 



£°(i,x,r-,p,7) := -TrKa°^ (f, x)7] + x)p + c°(t, x)r + /" (t, x) 
where the functions p, a, , b", c" and f" satisfy: 



|/x|oo + k|oo + SUp(|a°|i + |6"|i + |c"|l + |r|l) < 00. 

aeA 



Assumption HJB-|- The nonlinearity F satisfies HJB, and for any 6 > 0, there exists 
a finite set {aj}^l*^ such that for any a G A: 

inf |a--a"»|oo + |6°-6"1oo + |c"-c"^|oo + ir-r^|oo < S. 

l<i<Ms 



Remark 3.9. The assumption HJB+ is satisfied if ^ is a separable topological space and 

1.1 

(T°(-), c"(-) and /"(•) are continuous maps from A to C^' ; the space of bounded 

maps which are Lipschitz in x and ^"Holder in t. 

Theorem 3.10 (Rate of Convergence). Assume that the final condition g is bounded 
Lipschitz- continuous. Then, there is a constant C > such that: 
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(i) under Assumption HJB, we have v — < Ch^^'^ , 

(ii) under the stronger condition HJB+, we have —Ch^^^'^ < 



V — V 



The above bounds can be improved in some specific examples. See Subsection 3.4 for the 
hnear case where the rate of convergence is improved to y/h. 

We also observe that, in the PDE Finite Differences literature, the rate of convergence 
is usually stated in terms of the discretization in the space variable |Ax|. In our context 
of stochastic differential equation, notice that |Ax| is or the order of /i^/^. Therefore, the 
above upper and lower bounds on the rate of convergence corresponds to the classical rate 
|Aa;|^/^ and |Ax|^/^, respectively. 



3.2 Proof of the convergence result 



We now provide the proof Theorem 3.6 by building on Theorem 2.1 and Remark 2.1 of 
Barles and Souganidis [7j which requires the scheme to be consistent, monotone and stable. 
Moreover, since we are assuming the (weak) comparison for the equation, we also need to 
prove that our scheme produces a limit which satisfies the terminal condition in the usual 



sense, see Remark 3.2 



Throughout this section, all the conditions of Theorem 3.6 are in force. 



Lemma 3.11. Let ip be a smooth function with bounded derivatives. Then for all {t,x) S 
[0,r] X R<^: 

[c+ip]{t',x')-Th[c + ^]{t',x') 



lim 

(t' , x') — y (t, x) 
{h, c) -> (0, 0) 
t' + h <T 



h 



(£^^ + F(-,(^,Z)V^,Z)V)) {t,x) 



The proof is a straightforward application of Ito's formula, and is omitted. 



Lemma 3.12. Let ip^if) : [0,T] x M'^ — M be two Lip schitz functions. Then: 

^pKiP =^ Th[^]{t,x) <Th[ij]it,x) + ChE[{'ilj - ^){t + h,Xj;'')] for someC > 



where C depends only on constant K in (3.4). 



Proof. By Lemma |2.1| the operator can be written as: 



E 



V'CX^'") + HF [t,x,n'^{X]f)Hh{t,x)] 



Let f ■= ip — (f > where ip and ip are as in the statement of the lemma. Let Fj- denote 
the partial gradient with respect to r = {r,p,j). By the mean value Theorem: 



Thmt,x)-Th[ip]{t,x) 



E 
E 



/(X*'^') +hF^{9)-Vhf{X'^ 



t,X\ 



f{X];^){l + hFr{e)-Hn{t,x)) 



for some 9 = {t,x,r,p,^). By the definition of Hh{t,x): 

Thm-ThM = E [/(If) (1 + HFr + Fp.{a^r'Wh + h-^F^ ■ {a^)-\WhW^ - hl)a 
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where the dependence on 9 and x has been omitted for notational simphcity. Since < a 
by (3.4) of Assumption F, we have 1 — ■ > and therefore: 



= E 



nx'-){hF, + hF,.a^-''^' ' - ,T-i^.W^.^-i 



, +hF^-a' 
h 



/l2 



-a 



Let rup := max{— mi?,0}, where the function mp is defined in (3.4). Under Assumption 
F, we have K := \mp\cxi < oo, then 

F..a^-''^ + hF,.a^-''^a-^>-K 
h 



/l2 



one can write, 



Th[^]-Th[^] > E f{X]f){hFr-hK) 



> -C'hE 



t,X\ 



for some constant C > 0, where the last inequality follows from (3.4). 



□ 



The following observation will be used in the proof of Theorem 3.10 below 



Remark 3.13. The monotonicity result of the previous Lemma 3.12| is slightly different 
from that required in [7j. However, as it is observed in Remark 2.1 in [7], their convergence 
theorem holds under this approximate monotonicity. From the previous proof, we observe 
that if the function F satisfies the condition: 



Fr 



1 



(3.6) 



fp'F-Fp > 0, 
then, the standard monotonicity condition 

^ < il, ^ Th[ip\{t,x) < ThW\{t.x) (3.7) 

holds. Using the parabolic feature of the equation, we may introduce a new function 



u{t,x) := e^^^ ^^v{t,x) which solves a nonlinear PDE satisfying (3.6). Indeed, direct cal- 
culation shows that the PDE inherited by u is: 



C^u- F {■,u,Du,D'^u) =0, on [0,r)xM'^ 
u{T,x)=g{x), on M'', 



(3.8) 
(3.9) 



where F{t, x, r,p, 7) = e^'^'^~^^ F{t, x, e"^(^"*V, e"''^^"*)^, e"^('^"*)7) + Or. Then, it is easily 



seen that F satisfies the same conditions as F together with (3.6) for sufficiently large 6. 

■ M he two L°^- 



Lemma 3.14. Let ip,ip : [0,T] x M'^ 
exists a constant C > such that 



'bounded functions. Then there 



\Th[ip] - ThM^ < Iv? - V|oo(l + Ch) 



In particular, if g is E'^ — bounded, the family {v^)h defined in (2.4) is E^" — hounded, uni- 
formly in h. 
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Proof. Let f := (f — ip- Then, arguing as in the previous proof, 
where 



f{Xh) [l-a-^-F^ + h\Ah\^ + hFr - -rFjF-Fp 



bmce 1 - Tr[a-iF^] > 0, \Fr\oo < oo, and \F^F-Fp\ oo < oo by Assumption F, it follows 
that 

\Th[^]-Thm^ < \f\oo{l-a-' ■F^ + hE[\Ah\^] + Ch) 
But, = \F^F-Fp + a-^ ■ F^. Therefore, by Assumption F 

|T,M-T,[V;]L < \f\oo{l + \F^F-Fp + Ch^ <\f\Ul + Ch). 

To prove that the family {v^)h is bounded, we proceed by backward induction. By the 
assumption of the lemma v^{T, .) = g is L°°— bounded. We next fix some i < n and we 
assume that \v^{tj, .)|oo < Cj for every i + l<j<n — 1. Proceeding as in the proof of 
Lemma 3.12 with ip = v^{ti^i, .) and V' = 0, we see that 

v^iU,.) < h\F{t,x,0,0,0)\+Ci+i{l + Ch). 

oo 

Since F(t, x, 0,0,0) is bounded by Assumption F, it follows from the discrete Gronwall 
inequality that \v^{ti, .)\ao < Ce*^^ for some constant C independent of h. □ 



Remark 3.15. The approximate function defined by (2.4) is only defined on {ih\i = 
0, • • • ,N} X M"'. Our methodology requires to extend it to any t G [0,?"]. This can be 
achieved by any interpolation, as long as the regularity property of mentioned in Lemma 



3.16| below is preserved. For instance, on may simply use linear interpolation. 
Lemma 3.16. The function is Lipschitz in x, uniformly in h. 

Proof. We report the following calculation in the one-dimensional case d = 1 in order to 
simplify the presentation. 

1. For fixed t € [0, T — /i] , we argue as in the proof of Lemma 3.12 to see that for x, x' S 
with X > x': 

v^{t,x) -v^{t,x) = A + hB, (3.10) 
where, denoting S^^^i := D^v^{t + h, X\f) - D^v^{t + h, X]f^') for A: = 0, 1, 2: 

- F(t,x ,Vv^{t + h,X\l''') 



A := ¥.[5^^^]+h{F(t,x',Vv^{t + h,X];- 
= E [(1 + hFr)5^^^ + hFp5^^^ + ^ 



151 



F[t,x, Vv^it + h, X*'^) ) - F[t,x',Vv''{t + h, X*" 



\x — X \, 
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by Assumption F (i). By Lemma 2.1 we write for k = 1,2: 
E[5W] = E[6^^^H^it,x)+v\t + h,r/)(Hl:{t,x)-Hj^it,x') 



E[6'^°^Hj^{t,x)+Dv''{t + h,X'/] 



h 



k-l 

[a{t,x)-^ - a{t,x')-^'^cj{t,x')\. 



Then, dividing both sides of (3.10) by x — x' and taking hmsup, if follows from the above 
equalities that 

\v^(t,x) -v^{t,x')\ 
limsup 

|x-2;'|\0 [X — X ) 



< E 



lim sup 

|x-x'|\0 



[X — X' 



Wl - h 



v\t + h,X];-)-v\t + h,X\l^) ^ 



+Dv^{t + h,X];'') ( WhF-, 
\v^{t,x) - v^{t,x') 



'io-x{t,x) (T:^{t,x) 



^ a{t,xY ^ a{t,x 
2. Assume v^{t + /i, .) is Lipschitz with constant L^+z^. Then 

lim sup 



+ Ch. 



|a:-3;'|\0 ~ 



hN N'^ 

(1 + x)h + a^{t, x)VhN) [l + hFr + Fp^^—^ + F^ 



a{t,x) a{t,xY cr{t,xy 
+ Ch. 



Observe that 



p fx Fp \/F\ 



CJ 



'F^ a 



Since all terms on the right hand-side are bounded, under our assumptions, it follows that 
l-Pp^loo < CO (we emphasize that the geometric structure imposed in Assumption F (iii) 
provides this result in any dimension). Then: 

\v^(t,x) - v^(t,x')\ 



limsup 

|2:-x'|\0 [X — X) 



VhN 



Ar2 



(1 + ,At, + aAt, x)VhN) (1 + Fp— + F,-^ - 
+^NF, ~'^'''f\f ]+Ch]+ Ch. 
3. Let P be the probability measure equivalent to P defined by the density 



Z := 1 — a + aX"^ where 



a 



a{t,xy 



Then, 



\v^{t,x)-v^{t,x')\ 
limsup / ^ -<Lt+h\ ^ 

\x-x'WQ \X — X } 



1 + iix{t,x)h + a.j,{t,x)VhN') (l + Z^'^Fp 



hN 



a{t, x) 



a{t, x) 



+ Ch] + Ch. 



14 



By Cauchy-Schwartz inequality, we have 



\v''(t,x) - v'^(t,x')\ ^ 
limsup ' ^ ' ^ 7^^^<Lt+h 



\x-x'\\0 



X — X 



1 + fi^{t,x)h + a^{t,x)VhN)(l + Z ^Fp 



'hN 

' <T(t, x) 



a{t,xy 



+ Ch] +Ch 



By writing back the expectation in terms of probabiUty 

Z 



\v'^(t,x) - v'^(t,x')\ ^ 



\x-x'\\0 



X — X 



1 + ^^(t, x)h + a^{t, x)VhN) (l + Z-^Fp^^^ 

/ V cr(t, X) 



+Z-'VhNF^ 2a.(t x) 

By expanding the quadratic term inside the expectation, we observe that expectation of all 
the terms having ^/h, is zero. Therefore, 



\v''(t,x) - v'^(t,x')\ ^ 
limsup ' ^ '/ j^^^<Lt+h\ r 

|x-x'|\0 [X — X ) 



1 + ^^{t, x)h + a^{t, x)VhN) (l + Z-^Fp^^^ 

/ V O'it, X) 



+Z-^VhNF^ 



-2ax{t,x) 



^ (j{t,xf 
<Lt^h{{l + C'h)'^ +Ch] +Ch, 



+ Ch] +Ch 



which leads to 



lim sup 



\v^'{t,x) - v^{t,x') 



\x~x'\\0 [X - X') 

for some constants C, C > 0. 



□ 



Finally, we prove that the terminal condition is preserved by our scheme as the time step 
shrinks to zero. 

Lemma 3.17. For each x G M"' and tf^ = kh with k = 1, - ■ ■ ,n, we have; 

\v''{tk,x)-g{x)\ < C{T-tk)l 
Proof. 1. By the same argument as in the proof of Lemma |3. 14 we have: and for j > i: 

v\t,,X\:p = Et^ [v^{t,+,,Xl-^^) (1 - a, + a,Nf)_ 

where F^ := F{tj, Xj:^'^ , 0, 0, 0), aj, Fr , Fp are J^t^ — adapted random variables defined as in 



3.14 



at tj, and Nj 



the proof of Lemma 

Combine the above formula for j from i to n — 1, we see that 

n-l 



. —Wt- 

— ^-^^ — ^ has a standard Gaussian distribution. 



v\u,x)=E g{X'^''')P,,n +hEY,Fi+F}Et^[v\t,+^,X^;^;)] + Fl.Et^\Dv\t^+^,r^'^^^^^ 
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where Pi^k '■= \]!j=i ~ + '^i^'j) ^ ^ ^-^^ for all 1 < i < A; < n and Pi^i = 1. Obviously 
{Pi,kT'i' < A; < n} is a martingale for all i < n, a. property which will be used later. Since 



\F{-, •, 0, 0, 0)|oo < +00, and using Assumption F and Lemmas 3.16 and 3.14 

\vHti,x)-gix)\ < E 



giX'r^^'') - g{x)) P,^n +C{T-ti). 



(3.11) 



2. Let {(7e}e be the family of smooth functions obtained from g by convolution with a 
family of mollifiers {pe}, i.e. g^ = g * pe. Note that we have 



be - ^loo < Ce, \Dge\oo < \Dg\oo and \D'^g^\^ < e ^\Dg\^ 



(3.12) 



Then: 



E 



ff(X^-") - g{x) ) P, 



< E 
+ 



< Ce + 

< Ce + 



<7(X^"")-5e(^T''")^'i,n 

geiX^T^^) - ge{x)) Pi,n + {qs - g\ 

E 



E 



+ 



geiX'T^"") - ge{x)) Pi, 

E[Pi,„ r Dge{X';''')a{s)dWs\ 
Ju 



(3.13) 



where we denoted h{s) = b{tj, X^^.'^) and cr{s) = a{tj,X^^.'^) for tj < s < ij+i and d = a a 
We next estimate each term separately. 
2. a. First, since {Pi,k,i < k < n} is a martingale: 

pT n—l rt+i 

E[Pi,n Dge{Xl-'')a{s)dWs]\ = |^E[P,,„ / ' Dge{Xl-'')a{s)dWs] 



j=i 



< ^|e[P,-,+i / Dg,{Xl-^)a{s)dWs] 



n-1 



\E[P,,,ait,)Et^ / Dg.iX'/ndWs]] 

j=i -^^i 



Notice that 



E. 



P. 



■j+i 



Dg,{X'r)dWs 



E. 



E, 



W, 



■i+i 



Dg,{X'r)dWs 



■i+i 



2WsDg,{Xl-^)ds 
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Using Lemma 2.1 and (3.12), this provides: 



u 

n-1 



(3.14) 



j=i ■^ij 

n-1 

< Ce-^^h < C'{T-ti)e-\ 



(3.15) 



2.C. By (3.12) and the boundedness of b and a, we also estimate that: 

< C + Ce-\ 



I. 



(3.16) 



2.b. Plugging ( 3.15| ) and (3.16) into (3.13), we obtain: 



E 



QeiX'-'') - g,{x)) PiJ < C{T-ti) + C{T-ti)e-\ 



which by (3.11) provides: 

\v^{ti,x)-g{x)\ < Ce + C{T-ti)e'^ + C{T-ti). 



□ 



The required result follows from the choice e = ^/T — ti. 

Corollary 3.18. The function is 1/2-Hdlder continuous on t uniformly on h. 

Proof. The proof of ^-Holder continuity with respect to t could be easily provided by 
replacing g and v^(tk, •) in the assertion of Lemma respectively by v^{t, •) and v^{t' , ■) and 
consider the scheme from to time t' with time step equal to h. Therefore, we can write; 

\v^{t,x)-v^{t',x)\ < C{t'-t)^, 

where C could be chosen independent of t' for t' < T. 



□ 



3.3 Derivation of the rate of convergence 



The proof of Theorem 3.10 is based on Barles and Jakobsen [6], which uses switching 
systems approximation and the Krylov method of shaking coefficients 

3.3.1 Comparison result for the scheme 



Because F does not satisfy the standard monotonicity condition (3.7) of Barles and Sougani- 



dis [7], we need to introduce the nonlinearity F of Remark 3.13 so that F satisfies (3.6) 
Let u'^ be the familiy of functions defined by 



M^(r,.) =5 and u''(ti,x) =Th[u'']{ti,x), 
where for a function tp from [0, T] x M'^ to M with exponential growth: 

Thmt,^) :=E U(i + /i,^i'")l +hF{-,Vhij){t,x), 



(3.17) 
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and set 

v^{ti,x) := e~^^'^~^'^u''{ti,x), i = 0,...,n. (3.18) 

The following result shows that the difference — is of higher order, and thus reduces 
the error estimate problem to the analysis of the difference — v. 

Lemma 3.19. Under Assumption F, we have 

\imsuph~^\{v'^ -v^){t, .)\oo < oo. 

Proof. By definition of F, we directy calculate that: 

v''{t,x) = e'~^^{l + he)E[v''{t + h,X^^'')] + hF (t + h,x,VhV^{t,x] 



Since l + hO = e'^^ + OQi^), this shows that v^{t,x) = Th[v'^]{t, x) + 0{h'^). By lemma |3.14 
we conclude that: 

\{v^'-v''){t,-)\^ < {l + Ch)\{v^-v''){t + h,-)\^ + 0{h'), 

which shows by the Gronwall inequality that \{v^ — v^){t, •)|oo < 0{h) for all t < T — h. □ 



By Remark 3.13, the operator satisfies the standard monotonicity condition (3.7): 



^ < ^ ^ T^[ip] < Thii^]. (3.19) 

The key-ingredient for the derivation of the error estimate is the following comparison result 
for the scheme. 

Proposition 3.20. Let Assumption F holds true, and set /3 := |Fr|oo- Consider two 
arbitrary bounded functions (p and ip satisfying: 

h-^if -Thiif]) < 91 and h-^ -Th[ip]) > 92 (3.20) 
for some bounded functions gi and 92- Then, for every i = 0, • • • ,n: 

{^-m^,^) < e^(^-*'^l(</'-V')+(r,-)|oo + (T-/i)e''(^-*>)|(5i-<72)+|oo. (3.21) 

To prove this comparison result, we need the following strengthening of the monotonicity 
condition: 

Lemma 3.21. Let Assumption F hold true and let (3 := |-Fr|oo- Then, for every a,b € M+, 
and every bounded functions f < ip, the function 5{t) := e^^'^~^\a + b{T — t)) satisfies: 

Th[(p + S]{t,x) < Th[i^]{t,x) + 6{t)-hb, t<T-h,xeR'^. 
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Proof. Because 6 does not depend on x, we have V^lip + 6] = V^ip + 6{t + h)ei, where 
ei := (1, 0, 0). Then, it follows from the regularity of F that there exist some ^ such that: 

F{t + h, X, Vh [if + 5\ {t, x)) = F{t + h,x, Vhv{t, x)) + 6{t + h)Fr {t + h, x, + Vhf{t, x)) , 

and 

Th[f + 6]{t,x) = 6{t + h)+E[(f{t + h,xj;'')] + hF{t + h,x,Vhv{t,x)) 
+h6{t + h)Fr {t + h, X, ^ei + Vhfit, x)) 
= Th[(p] {t, x) + 5{t + h) {1 + hFr {t + h, X, Cei + Vhf{t, x)) } 
< Th[ip]{t,x) + {l + (3h)6{t + h). 



Since Th satisfies the standard monotonicity condition (3.19), this provides: 



Th[ip + 6]it,x) <Th[ij]it,x) + 6{t) + C{t), where C{t) ■= {1 + f3h) 5{t + h) - 5{t). 

It remains to prove that ((t) < —hb. From the smoothness of 5, we have 5{t + h) — 5{t) 
h6'{t) for some t £ [t,t + h). Then, since 5 is decreasing in t, we see that 

h-^C{t) = 6'{t}+ ^5{t + h) < 6'{t) + p5{t) < -6e^M, 

and the required estimate follows from the restriction 5 > 0. 



□ 



Proof of Proposition 3.20, We may refer directly to the similar result of f6]. However 
in our context, we give the following simpler proof. Observe that we may assume without 
loss of generality that 



ip{T,-) <7p(T,-) and gi < 92- 
Indeed, one can otherwise consider the function 

^ :=^ + e'^(^-*)(a + 6(r-t)) where a = | ((/? - 



(3.22) 



b = \{gi- 92) 



and /3 is the parameter defined in the previous Lemma 3.21 so that il:{T, •) > (/?(T, •) and 



by Lemma (3.21), ip{t,x) — ThWlii^x) > h{gi V (72)- Hence (3.22) holds true for 99 and tp 



We now prove the required result by induction. First ip{T, •) < t{j{T, •) by (3.22). We next 
assume that Lp(t + h,-) < ip{t + /i, •) for some t + h < T. Since Th satisfies the standard 



monotonicity condition (3.19), it follows from (3.22) that 



ThMt,x) < Thmt,x) 



On the other hand, under (3.22), the hypothesis of the lemma implies: 
ip{t,x)-ThMt,x) < ^{t,x)-ThMt,x). 

Then ip{t,-) < tl^it,-). 



□ 
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3.3.2 Proof of Theorem [37101 (i) 

Under the conditions of Assumption HJB on the coefficients, we may build a bounded 
subsolution ]/ of the nonhnear PDE, by the method of shaking the coefficients, which is 
Lipschitz in x, 1/2— Holder continuous in t, and approximates uniformly the solution v: 

V — e < v'^ < V. 

Let p{t,x) be a positive function supported in {{t,x) : t £ [0,1], |x| < 1} with unit 
mass, and define 



uf{t,x) := v'^ * where p^{t,x) := —j^p 



t X 



so that, from the convexity of the operator F, 



uf is a subsolution of (2.1), \uf — v\ < 2e. 
Moreover, since is Lipschitz in x, and 1/2— Holder continuous in t, 



(3.23) 



(3.24) 



vf is C^, and 



< Ce 



l-2/3o-|/3|i 



for any (/3o, /?) G N x N°' \ {0}, (3.25) 



where |/3|i := X^f^iA, and C > is some constant. As a consequence of the consistency 
result of Lemma |3. 11 above, we know that 

uf{t,x)-Th[w%t,x) 



nh[uf]{t,x) 



+ C^uf{t,x)+F{-,nf,D'uf,D^uf){t,x) 



h 



converges to as /i —t- 0. The next key-ingredient is to estimate the rate of convergence of 
TZhi'Uf] to zero: 



Lemma 3.22. For a family {(pe}o<e<i of smooth functions satisfying (3.25), we have: 

l^/iIVelloo — R{h, e) := C he^^ for some constant C > 0. 

The proof of this result is reported at the end of this section. From the previous estimate 
together with the subsolution property of w^, we see that uf < Th[uf] + Ch'^e-^. Then, it 
follows from Proposition 3.20 that 

uf -v'' < C\{yf - v^){T, .)|oo + Che-^ < C{e + he-^). (3.26) 

We now use (3.24) and ( |3.26| ) to conclude that 

V -v^ < V - nf + nf -v^ < C{e + he~^). 

Minimizing the right hand-side estimate over the choice of e > 0, this implies the upper 
bound on the error v — v^: 



V — V 



(3.27) 
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3.3.3 Proof of Theorem [sTTol (ii) 

The results of the previous section, together with the reinforced assumption HJB+, allow 
to apply the switching system method of Barles and Jakobsen [6] which provides the lower 
bound on the error: 

y-yh > -iniWe^^^ + R(h,e)} = -C'h^'^^, 

for some constants C, C" > 0. The required rate of convergence follows again from Lemma 
which states that the difference — is dominated by the above rate of convergence. 



3.19 



Proof of Lemma 



3.22 



Notice that the evolution of the Euler approximation X^jf" be- 
tween t and t + h\s driven by a constant drift ^{t, x) and a constant diffusion (j{t, x). Since 
D(pi; is bounded, it follows from Ito's formula that: 



^ipe{t + h, XI) - ^,{t, x)] - C^Mt, x) = 



t+h 



£^'-''ipe{u,X:)-C''ipe{t,x))du, 



where is the Dynkin operator associated to the Euler scheme: 

C^''''ip{t', x') = dtip{t', x') + fi{t, x)D(p{t',x') + ^Tr [a{t, x)D'^ip{t', x')] . 

Applying again Ito's formula, and using the fact that C^^ ^^Dip^ is bounded, leads to 



1 

h L" 



Eipe{t + h,X^) -ip,{t,x) -C^ipe{t,X 



1 



E 



t+h i-u 



C^''^C^'-''ip,{s,X^)dsdu. 



Using the boundedness of the coefficients /x and a, it follows from (3.25) that for e G (0, 1): 
E^,{t + h,X^)-ip,{t,x) 



C^ipe{t,x] 



h 



< Roih,e) :=C he~-\ 



Step 2 This implies that 

\Tlh\fe\{t.x)\ < 



Eips{t + h,Xj''')-ipe{t,x) 



h 



+ \F{X, ^eit, x),Dipeit, x),D\eit, x)) - F (•, P;,[99,] (t, x)) \ 
2 

< Ro{h, e) + CY^ |lE^Ve(t + h, X\f) - D^^e{t, x) 



k=0 



(3.28) 



by the Lipschitz continuity of the nonlinearity F. 
By a similar calculation as in Step 1, we see that: 

\ED'ipe{t + h,Xj;'') - Dipe{t,x)\ < Chs-^-\ i = 0,1,2, 



which, together with (3.28), provides the required result. 



□ 
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3.4 The rate of convergence in the hnear case 

In this subsection, we specialize the discussion to the Hnear one-dimensional case 

F(7) = C7, (3.29) 

for some c > 0. The multi-dimensional case d > 1 can be handled similarly. Assuming that 
g is bounded, the linear PDE (|2.1|)-(2.2 ) has a unique bounded solution 

M''. (3.30) 



v{t,x) = E[g{x + VlT2^WT^t)] for {t,x)£[0,T] 
We also observe that this solution v is C°° ([0,T) x M) with 



D^v{t,x) 



E 



g^^^ ix + Vl + Yc Wt^ 



, t < T, G M. 



(3.31) 



This shows in particular that v has bounded derivatives of any order, whenever the terminal 
data g is C°° and has bounded derivatives of any order. 

Of course, one can use the classical Monte Carlo estimate to produce an approximation 
of the function v of (3.30). The objective of this section is to analyze the error of the 
numerical scheme outlined in the previous sections. Namely: 



f ^(T, •) = 5, f '^(ti_i, x) = E U(ti, X + Wh) 



+ c/iE 



v^{ti,x + Wh)H!^] , i<n. (3.32) 



Here, a = 1 and /u = are used to write the above scheme. 



Proposition 3.23. Consider the linear F of {3.29), and assume that 15(2^+1) V is bounded 
for every k >0. Then 



limsup h ^^"^Iv^ — v\ 

h-s>0 



< oo. 



Proof. Since v has bounded first derivative with respect to x, it follows from Ito's formula 
that: 



v{t, x) 



E[v{t + h,x + Wh)] + cE 



[ Av{t + s,v + Ws)ds 
Jo 



Then, in view of Lemma 2.1 , the error u := v — v^ satisfies u{tn, Xt^) = and for i <n—V. 

u {ti ,XtJ = Ei [u {ti+i , Xt^^^ )]+chEi [An (t^+i , Xt^^^ ) ] 

+cEi [Aviih + s,Xih+s)-Av{{i + l)h,X(^i+^)h)]ds, (3.33) 
Jo 

where Ej := E[-|J^tJ is the expectation operator conditional on Tt^. 
Step 1 Set 



af :=E 



A^uit^Xt, 



:= E 



A'^v {ti^i + s, Xt^_,+s) - A^n (t„ XtJ 



ds. 
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and we introduce the matrices 

/ 1 -1 
1 -1 



A :-- 



\ 




V 



1 -1 

1 



/ 1 



B :-- 



\ 



\ 


1 

/ 



and we observe that (3.33) impHes that the vectors a'^ := (o^, . . . , a^)^ and := {bi, . . . , b'^)'^ 
satisfy Aa'' = chBa^^^ + cBb^ for all /c > 0, and therefore: 



1^ = chA-^Ba^^^ + cA-^Bb^ where ^"^ 



/l 1 
1 

Vo ... 



1 



1 



By direct calculation, we see that the powers {A ^B)^ are given by: 



j - 
k-1 



for all k > 1 and i,j = 1, . . . , n. 



In particular, because = 0, {A ^B)^ ^a^ = 0. Iterating (3.34), this provides: 



n-2 



ch{A-'B)a' +c{A-^B)b^ 



Y,c''^^h\A-^B)^+^b\ 



k=0 



and therefore: 



n-2 



fc=0 



Because of 



{A ^b)Ij - i{j>i+k} 



for all k > 1 and j = 1, . . . ,n , 



we can write (3.35): 



n—2 n 

u{0,x) = cY,{chf 

k=0 j=k+2 



By changing the order of the summations in the above we conclude that: 

n j-2 



u{0, x) 



j=2 k=0 



^-']b^-\ 
k J 



(3.34) 



(3.35) 



(3.36) 



Step 2 From our assumption that D^'^^^v is L°°— bounded for every A; > 0, it follows that 

cU 



ds 
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for some constant C. We then deduce from (3.36) that: 



n j-2 

(o,x)i < cch^/^ EE("^) 

j=2 k=0 



k 



So, 



J=2 



ch 



n-l 



< cVh. 



4 Probabilistic Numerical Scheme 



In order to implement the backward scheme (2.4), we stiU need to discuss the numerical 



computation of the conditional expectations involved in the definition of the operators 



in (2.5). In view of the Markov feature of the process X, these conditional expectations 



reduce to simple regressions. Motivated by the problem of American options in financial 
mathematics, various methods have been introduced in the literature for the numerical 
approximation of these regressions. We refer to |9] and [20] for a detailed discussion. 
The chief object of this section is to investigate the asymptotic properties of our suggested 



numerical method when the expectation operator E in (2.4) is replaced by some estimator 
corresponding to a sample size A^: 



E 



N 



^{t + h,Xl) +hF[-,VrM it,x), 



(4.1) 
(4.2) 



where 



Vhijit,x):=E^ ^{t + h,X]f)Hh{t,x) , Kh[^]:=U\\W^ + Cih) + C2K 



where 



Ci 



F F~ F \ -\- \F 



and C2 = |F(t,x,0,0, 



The above bounds are needed for technical reasons which were already observed in [9]. 
With these notations, the implementable numerical scheme is: 



v'^{t,x,u) = Tf^[v%]{t,x,io) 



(4.3) 



where is defined in (|4.1|)-(|4.2|), and the presence of u throughout this section emphasizes 



the dependence of our estimator on the underlying sample. 

Let TZh be the family of random variables R of the form 'ip(Wh)Hi(yVh) where ip is a 
function with iV'loo ^ b and Hi's are the Hermite polynomials: 



Ho{x) = 1, Hi{x) = X and H2{x) = x^ x - h Vx e R". 



Assumption E There exist constants Cb, A, > such that 
for every R G IZ^,, for some p > 1. 



¥.'^[R] -¥.[R\ 



< Cbh-^N- 
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Example 4.1. Consider the regression approximation based on the MaUiavin integration 
by parts as introduced in Lions and Reigner j27], Bouchard, Ekeland and Touzi [3J, and 
analyzed in the context of the simulation of backward stochastic differential equations by 
|U] and [13] • Then Assumption E is satisfied for every p > 1 with the constants A = and 
= ^, see [9]. 



Our next main result establishes conditions on the sample size A'" and the time step h 
which guarantee the convergence of towards v. 

Theorem 4.2. Let Assumptions E and F hold true, and assume that the fully nonlinear 
PDE (2.1) has comparison with growth q. Suppose in addition that 

lim/i^+^TV^ = oo. (4.4) 

Assume that the final condition g is bounded Lipschitz, and the coefficients fi and a are 
bounded. Then, for almost every lo: 

v^^{-,uj) — > V locally uniformly. 



where v is the unique viscosity solution o/(2.1). 



Proof. We adapt the argument of [7] to the present stochastic context. By Remark 3.13 



and Lemma 3.19, we may assume without loss of generality that the strict monotonicity 
( [3^ holds. 

By (4.2), we see that is uniformly bounded. So, we can define: 

v^{t,x):= liminf v^{t' ,x) and v*{t,x) := limsup v^{t' ,x'). (4.5) 

it' ,x') ^ (t,x) it' ,x') ^ {t,x) 

Our objective is to prove that v^, and v* are respectively viscosity superpersolution and 
subsolution of (2.1). By the comparison assumption, we shall then conclude that they 
are both equal to the unique viscosity solution of the problem whose existence is given by 
Theorem 3.6 In particular, they are both deterministic functions. 

We shall only report the proof of the supersolution property, the subsolution property 
follows from the same type of argument. 

In order to prove that v^, is a supersolution of (2.1), we consider {tQ,XQ) G [0,T) x M" 
together with a test function (/? G ([0,T) x M"), so that 

= min{{)* = {v* - ip){to,xo). 

By classical manipulations, we can find a sequence (t„, x„, /i„) — )• {t^, xq, 0) so that v^"^ (t„, Xn) 
v^,{tQ,xo) and 

{v^" - if){tn, Xn) = min{{)''" - 99} =: Cn ^ 0. 
Then, ■O'^" > 9? + Cn, and it follows from the monotonicity of the operator that: 

T,J^'^"]>T,J(^ + C„]. 
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By the definition of w'^" in (4.3), this provides: 

v'^- {t,x)> [ip + Cn] {t,x)- (T,„ - t,„ ) it] (t, x) , 

where, for ease of notations, the dependence on has been dropped. Because v'^" {tn, x„) = 
ip{tn, Xn) + Cn, the last inequality gives: 

^{tn,Xn) + Cn- ThA'P + Cn]{tn,Xn) + KRn > 0, Rn ■= K^{Th„ - T hj[v''"]{tn, Xn) ■ 

We claim that 

Rn — > P — a.s. along some subsequence. (4-6) 

Then, after passing to the subsequence, dividing both sides by hn, and sending n — )• oo, it 
follows from Lemma 13.111 that : 

-C^ip - F (•, ip, Dip, > 0, 

which is the required supersolution property. 



It remains to show (4.6). We start by bounding R^ with respect to the error of estimation 
of conditional expectation. By Lemma 3.14, |T/i^ |oo < K^^^ and so by (4.2), we can 
write: 



< 



Th^-Th„)[v''-]itn,Xn) 



By the Lipschitz-continuity of F, we have: 

(T;,„ - t,,„) Xn)\ < C {So + K£l + hn£2) 



(4.7) 



where: 



|(E - ^)[v^-{tn + K,Xll)H'l-{tn,Xn)]\ 



Th^-ThA[v''-]{tn,Xn) 



< c 



where i?^ = v^" {tn + hn, Xn + a{x)Wh)Hi{Wh), i = 1, 2, 3 and Hi is Hermite polynomial of 
de gree i. This leads the following estimate for the error Rn- 



\Rn\ < 



c 



(E-]E)[i?0] + {E-E)[Ri 



+ h 



-1 



[E-t)[Rl] 



(4.8) 



Because R\ G IZb with bound obtained in Lemma 3.14 by Assumption E we have,: 



SO by (4.4) we have ||-Rn||p — > which implies (4.6). 



□ 



We finally discuss the choice of the sample size so as to keep the same rate for the error 
bound. 
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Theorem 4.3. Let the nonlinearity F be as in Assumption HJB, and consider a regression 
operator satisfying Assumption E. Let the sample size be such that 



Urn To Ni: > 0. 

h->0 



(4.9) 



Then, for any bounded Lip schitz final condition g, we have the following U' —bounds on the 
rate of convergence: 



\v — V 



< Ch'/''. 



Proof. By Remark 3.13 and Lemma 3.19 we may assume without loss of generality that 
the strict monotonicity (3.6) holds true. 



We proceed as in the proof of Theorem 3.10 to see that 

v-v^ < V -v'' + Vh- = e + R{h,e) + Vh- 
Since satisfies (4.3), 



h 



T^hlA] > -Rh[v''] where Rh[ip] := 



h 



where, in the present context, Rh[v^] is a non-zero stochastic term. By Proposition 3.20 
it follows from the last inequality that: 



v-v^ <C{e + R{h, e) + Rhiv""] ) , 



where the constant C > depends only on the Lipschitz coefficient of F, /3 in Lemma 3.21 
and the constant in Lemma 13.221 

Similarly, we follow the line of argument of the proof of Theorem 3.10 to show that a 
lower bound holds true, and therefore: 



\v — V 



<c(^e^/^ + R{h,e) + Rh[v^]) , 



We now use (4.9) and proceed as in the last part of the proof of Theorem 4.2 to deduce 
from (4.8) and Assumption F that 

\\Rh[A\\p < Ch'l^"^. 
With this choice of the sample size A^, the above error estimate reduces to 

\\v^-v\\p < C [e^/^ + R{h,e) + h^/^^' 

and the additional term h^^^^ does not affect the minimization with respect to e. □ 

Example 4.4. Let us illustrate the convergence results of this section in the context of the 
Malliavin integration by parts regression method of |27| and [9J where A = ^ and = ^ 
for every p > 1. So, for the convergence result we need to choose Nh of the order of h~"'' 
with ao > 5 +4p. For the L^-rate of convergence result, we need to choose of the order 
of with ai > ^ + ^. 
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5 Numerical Results 



In this section, we provide an application of the Monte Carlo-finite differences scheme sug- 
gested in this paper in the context of two different types of problems. We first consider 
the classical mean curvature flow equation as the simplest front propagation example. We 
test our backward probabilistic scheme on the example where the initial data is given by a 
sphere, for which an easy explicit solution is available. A more interesting geometric exam- 
ple in space dimensions 2 is also considered. We next consider the Hamilton-Jacobi-Bellman 
equation characterizing the classical optimal investment problem in financial mathematics. 
Here, we again test our scheme in dimension two where an explicit solution is available, and 
we consider more involved examples in space dimension 5, in addition to the time variable. 

In all examples considered in this section the operator F{t,x,r,p,j) does not depend on 
the r— variable. We shall then drop this variable from our notations, and we simply write 
the scheme as: 

v^{T,.):=9 and 

v^{ti,x):=nv''{ti+i,Xl)] + hF{U,x,VhV^{%x)) ^ ' ^ 

where 

We recall from Remark 



and Vl and Vl are defined in Lemma 



2.1 



2.2 



that: 



(5.2) 



Vly^ip{ti,x) = E 

= E 

The second representation is the one reported in jT2| where the present backward probabilis- 
tic scheme was first introduced. These two representations induce two different numerical 
schemes because once the expectation operator E is replaced by an approximation , 
equality does not hold anymore in the latter equation for finite N . In our numerical ex- 
amples below, we provide results for both methods. The numerical schemes based on the 
first (resp. second) representation will be referred to as scheme 1 (resp. 2). An important 
outcome of our numerical experiments is that scheme 2 turns out to have a significantly 
better performance than scheme 1. 

Remark 5.1. The second scheme needs some final condition for 'Dj^ip{T, Xl^~^'^). Since g 
is smooth in all our examples, we set this final condition to Vgf. Since the second scheme 
turns out to have a better performnace, we may also use the final condition for Z suggested 
by the first scheme. 

We finally discuss the choice of the regression estimator in our implemented examples. 
Two methods have been used: 

• The first method is the basis projection a la Longstaff and Schwartz [28], as devel- 
oped in j20j . We use regression functions with localized support : on each support 
the regression functions are chosen linear and the size of the support is adaptative 
according to the Monte Carlo distribution of the underlying process. 
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• The second method is based on the Mahiavin integration by parts formula as suggested 
in [27j and further developed in |3J. In particular, the optimal exponential localization 
function (t)^{y) = exp{—7]'^y) in each direction k is chosen as follows. The optimal 
parameter r]k is provided in and should be chosen for each conditional expectation 
depending on k. Our numerical experiments however revealed that such optimal 
parameters do not provide sufficiently good performance, and more accurate results 
are obtained by choosing r]k = 5/\/A^ for all values of k. 

5.1 Mean curvature flow problem 

The mean curvature flow equation describes the motion of a surface where each point moves 
along the inward normal direction with speed proportional to the mean curvature at that 
point. This geometric problem can be characterized as the zero-level set S{t) := {x G M*^ : 
v{t,x) = 0} of a function v{t,x) depending on time and space satisfying the geometric 
partial differential equation: 

^ Dv ■ D'^vDv 1 / \ / \ / N 

vt — Av -\ 1 p5 = and v[u,x) = g[x) (5-3) 

\Dv\'' 

and g :R'^ — > M is a bounded Lipschitz-continuous function. We refer to for more 
details on the mean curvature problem and the corresponding stochastic representation. 

To model the motion of a sphere in M'^ with radius 2R > 0, we take g{x) := AR^ — |xp so 
that g is positive inside the sphere and negative outside. We first solve the sphere problem 
in dimension 3. In this case, it is well-known that the surface S{t) is a sphere with a 



radius R{t) = 1\/R? — t for t G (0,i?^). Reversing time, we rewrite (5.3) for t € (0, T) with 
T = R^: 

-vt-^a'^Av + F{x,Dv,D'^v) = and viT,x) = g{x), (5.4) 

where 

F(x,z,7) := 7Q^'-l)+^- 

We implement our Monte Carlo-finite differences scheme to provide an approximation v'^ of 
the function v. As mentioned before, we implement four methods: Malliavin integration by 
parts-based or basis projection-based regression, and scheme 1 or 2 for the representation 
of the Hessian. 

Given the approximation v^, we deduce an approximation of the surface S^{t) := {x G 
: v^{t,x) = 0)} by using a dichotomic gradient descent method using the estimation 
of the gradient V^v estimated along the resolution. The dichotomy is stopped when the 
solution is localized within 0.01 accuracy. 

Remark 5.2. Of course the use of the gradient is not necessary in the present context 
where we know that S{t) is a sphere at any time t £ [0,T). The algorithm described above 
is designed to handle any type of geometry. 
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Remark 5.3. In our numerical experiments, the nonlinearity F is truncated so that it is 
bounded by an arbitrary value taken equal to 200. 

Our numerical results show that Malliavin and basis projection methods give similar 
results. However, for a given number of sample paths, the basis projection method of |20j 
are slightly more accurate. Therefore, all results reported for this example correspond to 
the basis projection method. 

Figure [T] provides results obtained with one million particles and 10 x 10 x 10 mesh with 
a time step equal to 0.0125. The diffusion coefficient a is taken to be either 1 or 1.8. 
We observe that results are better with a = 1. We also observe that the error increases 
near time 0.25 corresponding to an acceleration of the dynamics of the phenomenon, and 
suggesting that a thinner time step should be used at the end of simulation. 



Hean Curvature FIdh in 3D for a sphere 




' ' ' ' ' ' 

e D,B5 0,1 0,15 0,2 0,25 

tine 



Figure 1: Solution of the mean curvature flow for the sphere problem 

Figure [2] plots the differ ence between our calculation and the reference for scheme 1 and 
volatility 1 and 1.8 for varying time step. The corresponding results with scheme 2 are 
reported in figure [3j We notice that some points at time T = 0.25 are missing due to a non 
convergence of the gradient method for a diffusion a = 1.8. We observe that results for 
scheme 2 are slightly better than results for scheme 1. With o" = 1, it takes 150 seconds on 
a Nehalem intel processor 2.9 GHz to obtain the result at time t = 0.15 with the regression 
method, while it takes 1500 seconds with the Malliavin method (notice that the dichotomy 
used with the gradient method is a very inefficient method). 

We finally report in Figure [4] some numerical results for the mean curvature flow problem 
in dimension 2 with a more interesting geometry: the initial surface (i.e. the zero-level 
set for v) consists of two disks with unit radius, with centers positioned at -1.5 and 1.5 
and connected by a stripe of unit width. We give the resulting deformation with scheme 2 
for a diffusion a = 1, a time step h = 0.0125, and one million particles. Once again, the 
Malliavin integration by parts based regression method and the basis projection method 
with 10 X 10 meshes produce similar results. We used 1024 points to describe the surface. 
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Figure 2: Mean curvature flow problem for difl'erent time step and diffusion: scheme 1 




Figure 3: Mean curvature flow problem for difl'erent time step and diff'usions: scheme 2 



One advantage of this method is the total parallelization that can be performed to solve 
the problem for different points on the surface : for the results given parallelization by 
Message Passing (MPI) was achieved. 

5.2 Continuous-time portfolio optimization 

We next report an application to the continuous-time portfolio optimization problem in 
financial mathematics. Let {St, t £ [0,T]} he an ltd process modeling the price evolution of 
n financial securities. The investor chooses an adapted process {Ot,t £ [0, T]} with values 
in M", where 61 is the amount invested in the i— th security held at time t. In addition, the 
investor has access to a non-risky security (bank account) where the remaining part of his 
wealth is invested. The non-risky asset is defined by an adapted interest rates process 
{rt,t £ [0,T]}, i.e. dS^ = S^rtdt, t G [0,1]. Then, the dynamics of the wealth process is 
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Figure 4: Mean curvature flow problem in 2D 



described by: 



dx! = 9,-^ + ix!-e,.i)^-^ = e,.^ + ix!-e,.i)ndt, 

where 1 = (1, • • • , 1) G M'^. Let A be the collection of all adapted processes 6 with values 
in M.'^, which are integrable with respect to S and such that the process is uniformly 
bounded from below. Given an absolute risk aversion coefficient rj > 0, the portfolio 
optimization problem is defined by: 

vq := supE — exp(— r/X|>) . (5-5) 

Under fairly general conditions, this linear stochastic control problem can be characterized 
as the unique viscosity solution of the corresponding HJB equation. The main purpose 
of this subsection is to implement our Monte Carlo-finite differences scheme to derive an 
approximation of the solution of the fully nonlinear HJB equation in non-trivial situations 
where the state has a few dimensions. We shall first start by a two-dimensional example 
where an explicit solution of the problem is available. Then, we will present some results 
in a five dimensional situation. 

5.2.1 A two dimensional problem 

Let d = 1, rt = for all t € [0, 1], and assume that the security price process is defined by 
the Heston model |2T]: 



dSt = fiStdt + y^tStdWt^^'^ 

dYt = k{m - Yt)dt + c^/Vt (pdW^^^^ + y^l - p'^dWf 
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where W = {W'^^\W^^'^) is a Brownian motion in M^. In this context, it is easily seen that 



the portfoho optimization problem (5.5) does not depend on the state variable s. Given 
an initial state at the time origin t given by {Xt^Yt) = {x,y), the value function v{t,x,y) 
solves the HJB equation: 

v{T, X, y) = -e'"^"^ and = -vt - k{m - y)vy - \c^yvyy - sup i^O^yVxx + OifJ^v^ + pcyv^y) 

em 

= -Vt - k{m - y)vy - \c^yVyy + 



^yvxx 

A quasi explicit solution of this problem was provided by Zariphopoulou 



1 r 1^ 



exp(--y^ ^ds 



T ,,2 



v{t,x,y) = -e ''^ 

where the process Y is defined by 

Yt = y and dY = {k{m - Yt) - ficp)dt + c\/ YtdWt- 



(5.6) 
(5.7) 



In order to implement our Monte Carlo- finite differences scheme, we re- write (5.6) as: 



-Vt - k{m - y)vy - ^c^yvyy - ^a'^Vxx + F (y, Dv, D^v) = 0, v{T, x, y) = -e"''^, (5.8) 
where cr > and the nonlinearity F : M x x §2 is given by: 

F[y,z,-f) = -o-7iiH . 

2 22/711 

Notice that the nonlinearity F does not to satisfy Assumption F, we consider the truncated 
nonlinearity: 

Fe,M{y,z,l) ■■= ^0-^711- sup (IrO^iyV e)-fii + e{pzi + pc{y\/ e)ji2] , 

^ e<e<M \^ J 

for some e, n > jointly chosen with a so that Assumption F holds true. Under this form, 
the forward two-dimensional diffusion is defined by: 



dXf^=odWl'', and dX\'' = k(m - X\^^)dt ^ X\'' dWr . (5.9) 

In order to guarantee the non-negativity of the discrete-time approximation of the process 
X^^^ we use the implicit Milstein scheme [22] : 



_ Xr_^, + A:mAt + cVx;,!^iWAt+|c^A(e^-l) 

where (^n)n>i is a sequence of independent random variable with distribution N(0, 1). 

Our numerical results correspond to the following values of the parameter: p = 0.15, 
c = 0.2, k = 0.1, m = 0.3, Yq = m, p = 0. The initial value of the portfolio is xq = 1, 
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the maturity T is taken equal to one year. With this parameters, the value function is 
computed from the quasi-explicit formula (5.7) to be vq = —0.3534. 

We also choose M = 40 for the truncation of the nonlinearity. This choice turned out to 
be critical as an initial choice of M = 10 produced an important bias in the results. 

The two schemes have been tested with the Malliavin and basis projection methods. The 
latter was applied with 40 x 10 basis functions. We provide numerical results correspond- 
ing to 2 millions particles. Our numerical results show that the Malliavin and the basis 
projection methods produce very similar results, and achieve a good accuracy: with 2 mil- 
lions particles, we calculate the variance of our estimates by performing 100 independent 
calculations: 

• the results of the Malliavin method exhibit a standard deviation smaller than 0.005 
for scheme one (except for a step equal to 0.025 and a volatility equal to 1.2 where 
standard deviation jumped to 0.038), 0.002 for scheme two with a computing time of 
378 seconds for 40 time steps, 

• the results of the basis projection method exhibit a standard deviation smaller than 
0.002 for scheme 1 and 0.0009 for scheme two with a computing time of 114 seconds 
for 40 time steps. 

Figure [5] provides the plots of the errors obtained by the integration by parts-based 
regression with Schemes one and two. All solutions have been calculated as the average 
of 100 calculations. We first observe that for a small diffusion coefficient a = 0.2, the 
numerical performance of the algorithm is very poor: surprisingly, the error increases as 
the time step shrinks to zero and the method seems to be biased. This numerical result 
hints that the requirement that the diffusion should dominate the nonlinearity in Theorem 



3.6, might be a sharp condition. We also observe that scheme one has a persistent bias 



Error for schene t 



< f inancial probli 



Error for schene tuo> financial problei 




Figure 5: Difference between calculation and reference for scheme one and two 



even for a very small time step, while scheme two exhibits a better convergence towards 
the solution. 
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5.2.2 A five dimensional example 

We now let n = 2, and we assume that the interest rate process is defined by the Ornstein- 
Uhlenbeck process: 



dn 



K{b - rt)dt + C,dW^ 



(0) 



While the price process of the second security is defined by a Heston model, the first 
security's price process is defined by a CEV-SV models, see e.g. [29j for a presentation of 
these models and their simulation: 



ds: 



lx,S?dt + a,JYf^sf'dwl:'^^\ /32 = 1, 



dY^^ = h {mi - y/^^) dt + Ci\J yI'Uw}''^^ 



where {W^°\w(^''\w(^'^\ W^^'^K VF(2.2)) 

is a Brownian motion in M^, and for simplicity 
we considered a zero-correlation between the security price process and its volatility process. 

Since /32 = 1, the value function of the portfolio optimization problem (5.5) does not 
depend on the s*^^^— variable. Given an initial state {Xt, rt, sI:^\y^^\ y/^^) = (x, r, si, yi, 7/2) 
at the time origin t, the value function v (t, x, r, si,yi, 2/2) satisfies the HJB equation: 

-Vt - (L"" + + L^^)v - rxvx 

- sup <9i- rl)vx + Oiajyisl'^'^'^v^si + ^(^I'^^yiSi'^'"^ + Olaly2)v^ 

01,92 I ^ 

Y 







+ 



Vt - (L"" + + L-^ )v - rxvx 
{{fil - r)v^ + cTjyisl^'-\^si? , ((^2 



2afyis 



2ft -2 
1 ^ 



+ 



i2 - r Vx 



2o■|y2^'^ 



(5.11) 



where 



1 ^ 1 



ql 

and L v 



Aiisi^si - -OiSiyiVs^si- 



In order to implement our Monte Carlo-finite differences scheme, we re-write (|5.11) as: 

(5.12) 



-vt - {17- + 1/ + l.^')v - \a\xx + F {{x, r, si,yi,y2), Dv, D%) = 0, 
v{T,x,r,si,yi,y2) = 



-r)x 



where cr > 0, and the nonlinearity F : 



F{u,z,'y) 



1 



X §2 is given by: 

2/3i-l 



-a 711 - X1X2Z1 + 



{{^ii - X2)zi + a(x4X.p 71,3)2 {{1^2 - X2)zi) 



+ 



2" ' 2ajx4xf'-'^ju ' 2^2x5711 

where u = {xi, ■ ■ ■ ,^5). We next consider the truncated nonlinearity: 

Fs^m{u,z,j) := }:a^jii- X1X2Z1+ sup \ {6 ■ {fi - rl)zi + Oiafix^V e){x3 V e)'^'^'~'^ 
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+l{(^Wi{x3 V e)(x4 V ef"'-^ + Qlalixr, V e))7ii}, 
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where e,M > are jointly chosen with a so that Assumption F holds true. Under this 
form, the forward two-dimensional diffusion is defined by: 

dxf ^ = adwf ^ , dX^ = K{b - X[^^)dt + Qdwi'^^ , 



dxf ) = /.iXf )dt + a^^xi'^^'dwi'''\ dxf ) = Mm, - Xi'^)dt + c,^/xi^ dwi'''\ 
dxP = k2{m2 - Xf'^)dt + C2^X^Uwl'^''^\ 

(5.13) 

(2) 

The component X]- is simulated according to the exact discretization: 



where (^ji)n>i is a sequence of independent random variable with distribution N(0, 1). The 
following scheme for the price of the asset guarantees non- negativity (see [I] ) : 



f 1 2) 21 f 1 2) 

where AVF„ ' := Wn ' — W^nJi • We take the following parameters fii = 0.10, cJi = 0.3, 
(3i = 0.5 for the first asset, ki = 0.1, mi = 1., ci = 0.1 for the diffusion process of the 
first asset. The second asset is defined by the same parameters as in the two dimensional 

(2) 

example: fi2 = 0.15, C2 = 0.2, m = 0.3 and Yq = m. As for the interest rate model we 
take b = 0.07, X^^^ = b, ( = 0.3. 

The initial values of the portfolio the assets prices are all set to 1. For this test case we 
first use the basis projection regression method with 4x4x4x4x10 meshes and three 
millions particles which, for example, takes 520 seconds for 20 time steps. Figure |6] contains 
the plot of the solution obtained by scheme 2, with different time steps. We only provide 
results for the implementation of scheme 1 with a coarse time step, because the method 
was diverging with a thinner time step. We observe that there is still a difference for very 
thin time step with the three considered values of the diffusion. This seems to indicate that 
more particles and more meshes are needed. While doing many calculation we observed 
that for the thinner time step mesh, the solution sometimes diverges. We therefore report 
the results corresponding to thirty millions particles with 4x4x4x4x40 meshes. First 
we notice that with this discretization all results are converging as time step goes to zero: 
the exact solution seems to be very closed to —0.258. During our experiments with thirty 
millions particles, the scheme was always converging with a very low variance on the results. 
A single calculation takes 5100 seconds with 20 time steps. 

Remark 5.4. With thirty millions particles, the memory needed forced us to use 64-bit 
processors with more than four gigabytes of memory. 



5.2.3 Conclusion on numerical results 

The Monte Carlo-Finite Differences algorithm has been implemented with both schemes 
suggested by (5.2), using the basis projection and Malliavin regression methods. Our 
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numerical experiments reveal that the second scheme performs better both in term of results 
and time of calculation for a given number of particles, independently of the regeression 
method. 

We also provided numerical results for different choices of the diffusion parameter in the 
Monte Carlo step. We observed that small diffusion coefficients lead to poor results, which 
hints that the condition that the diffusion must dominate the nonlinearity in Assumption 
F (iii) may be sharp. On the other hand, we also observed that large diffusions require 
a high refinement of the meshes meshes, and large number of particles, leading to a high 
computational time. 

Finally, let us notice that a reasonable choice of the diffusion could be time and state 
dependent, as in the classical importance sampling method. We have not tried any ex- 
periment in this direction, and we hope to have some theoretical results on how to choose 
optimally the drift and the diffusion coefficient of the Monte Carlo step. 
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